x_vec(year) = x;
m_vec(year) = m;
omega_vec(year) = omega;
%relative_product_wage_vec(year) = real_product_wage.US/real_product_wage.OECD;
relative_product_wage_vec(year) = real_product_wage.US/real_product_wage.OECD;
std_relative_qualities_vec(year) = std(qualities_vec.US./qualities_vec.OECD);
steals_US(year) = sum(maybe_stolen.US);
steals_OECD(year) = sum(maybe_stolen.OECD);
    
for cc=1:length(place)
    wage_level_vec.(place{cc}) = [wage_level_vec.(place{cc}); wage_level.(place{cc})];
    real_consumption_vec.(place{cc}) = [real_consumption_vec.(place{cc}); real_consumption.(place{cc}).(year_index)];
    
    % DETERMINE WHICH VARIETIES ARE OWNED BY ENTRANTS AND ASSIGN THEM NEW FIRM NUMBERS
    actual_entrants.(place{cc}).(year_index) = entrants.(place{cc}) & produced.(place{cc}).(year_index);
    varown_vec.(place{cc}).(year_index)(actual_entrants.(place{cc}).(year_index)) = (max(firm_list.(place{cc}))+1:max(firm_list.(place{cc}))+sum(actual_entrants.(place{cc}).(year_index)));
    actual_entrants_num_index.(place{cc}).(year_index) = (max(firm_list.(place{cc}))+1:max(firm_list.(place{cc}))+sum(actual_entrants.(place{cc}).(year_index)))';
    firm_list.(place{cc}) = [firm_list.(place{cc}); (max(firm_list.(place{cc}))+1:max(firm_list.(place{cc}))+sum(actual_entrants.(place{cc}).(year_index)))'];
 
    % ENTRANTS THAT DO NOT ACTUALLY START PRODUCING (not relevant in base
    % case)
    placeholder_entrants.(place{cc}) = entrants.(place{cc}) & ~actual_entrants.(place{cc}).(year_index);
    varown_vec.(place{cc}).(year_index)(placeholder_entrants.(place{cc}))=NaN;
    
    % MAKE SURE DORMANT VARIETIES THAT ARE IMPORTED DO NOT HAVE FIRM OWNERSHIP
    varown_vec_dormant.(place{cc})(produced.(place{cc}).(year_index)) = NaN;
    varown_vec_dormant.(place{cc})(~produced.(place{cc}).(year_index)) = varown_vec.(place{cc}).(year_index)(~produced.(place{cc}).(year_index));
    varown_vec.(place{cc}).(year_index)(~produced.(place{cc}).(year_index)) = NaN;
    dormant_still_alive.(place{cc}) = ismember(varown_vec_dormant.(place{cc}),varown_vec.(place{cc}).(year_index));
        
    % FIND VARIETIES THAT ARE PRODUCED BUT UNOWNED (BECAUSE THEY WERE NOT
    % INNOVATED UPON) AND ASSIGN THESE TO NEW ENTRANTS OR THE STILL LIVING
    % INCUMBENT WITH OTHER NON-DORMANT VARIETIES
    gain_by_relative_wages_short.(place{cc}) = isnan(varown_vec.(place{cc}).(year_index)(produced.(place{cc}).(year_index)));
    gain_by_relative_wages.(place{cc}) = false(size(varown_vec.(place{cc}).(year_index)));
    gain_by_relative_wages.(place{cc})(produced.(place{cc}).(year_index))=gain_by_relative_wages_short.(place{cc});
    varown_vec.(place{cc}).(year_index)(gain_by_relative_wages.(place{cc}) & dormant_still_alive.(place{cc})) = ...
                    varown_vec_dormant.(place{cc})(gain_by_relative_wages.(place{cc}) & dormant_still_alive.(place{cc}));
    varown_vec.(place{cc}).(year_index)(gain_by_relative_wages.(place{cc}) & ~dormant_still_alive.(place{cc})) = ...
                    (max(firm_list.(place{cc}))+1:max(firm_list.(place{cc}))+sum(gain_by_relative_wages.(place{cc})&~dormant_still_alive.(place{cc})))';
    
    % UPDATE FIRM LIST BY ADDING NEW ENTRANTS TO THE VECTOR
    firm_list.(place{cc})=[firm_list.(place{cc});(max(firm_list.(place{cc}))+1:max(firm_list.(place{cc}))+sum(gain_by_relative_wages.(place{cc})))'];    
    total_entrants.(place{cc}) = gain_by_relative_wages.(place{cc}) & ~dormant_still_alive.(place{cc}) | actual_entrants.(place{cc}).(year_index);    

    % TRACK THE YEAR IN WHICH EACH FIRM EXITED
        firm_exit.(place{cc})=[firm_exit.(place{cc});nan(sum(total_entrants.(place{cc})),1)];
        new_exit.(place{cc}) = ~ismember(firm_list.(place{cc}),varown_vec.(place{cc}).(year_index)(produced.(place{cc}).(year_index))) & isnan(firm_exit.(place{cc}));
        new_exit_num_index.(place{cc}).(year_index) = firm_list.(place{cc})(new_exit.(place{cc}));
        firm_exit.(place{cc})(new_exit.(place{cc})) = year+1;
               
    % TRACK FIRM BIRTH AND CURRENTLY LIVING FIRMS
        firm_birthday.(place{cc}) = [firm_birthday.(place{cc}); [year+zeros(sum(total_entrants.(place{cc})),1)]];
        currently_alive.(place{cc}).(year_index)=isnan(firm_exit.(place{cc}));
    
    % VARIETY OWNERSHIP OF EXPORTS
    varown_exports.(place{cc}).(year_index) = zeros(length(varown_vec.(place{cc}).(year_index)),1);
    varown_exports.(place{cc}).(year_index)(exports.(place{cc}).(year_index)) = varown_vec.(place{cc}).(year_index)(exports.(place{cc}).(year_index));
    
    % TRACK THE NUMBER OF VARIETIES AND EXPORTS OWNED BY EACH FIRM
    firm_produced.(place{cc}) = accumarray(varown_vec.(place{cc}).(year_index)(produced.(place{cc}).(year_index)),countvars.(place{cc})(produced.(place{cc}).(year_index)));
    firm_number_of_exports.(place{cc}) = accumarray(varown_vec.(place{cc}).(year_index)(exports.(place{cc}).(year_index)),countvars.(place{cc})(exports.(place{cc}).(year_index)));
    stacked_ownership=[varown_vec.(place{cc}).(year_index)(produced.(place{cc}).(year_index));varown_vec.(place{cc}).(year_index)(exports.(place{cc}).(year_index))];
    stacked_countvars=ones(size(stacked_ownership));
    firm_size.(place{cc}) = accumarray(stacked_ownership,stacked_countvars);
    
    % CALCULATE REVENUES AND OUTPUT STATS FOR EACH FIRM
    stacked_revenue = [revenue.(place{cc}).(year_index)(produced.(place{cc}).(year_index)); ...
                        revenue.(place{3-cc}).(year_index)(exports.(place{cc}).(year_index))/import_tariff.(place{3-cc})];                               
    stacked_output = [output.(place{cc}).(year_index)(produced.(place{cc}).(year_index)) ; ...
                                   output.(place{3-cc}).(year_index)(exports.(place{cc}).(year_index))];                          
    firm_revenue.(place{cc}).(year_index) = accumarray(stacked_ownership,stacked_revenue);
    firm_output.(place{cc}).(year_index) = accumarray(stacked_ownership,stacked_output);
    std_dev_firm_output.(place{cc})(end+1) = nanstd(log(firm_output.(place{cc}).(year_index)(firm_output.(place{cc}).(year_index)~=0)));
    firm_export_output.(place{cc}).(year_index) = accumarray(varown_vec.(place{cc}).(year_index)(exports.(place{cc}).(year_index)),...
                            output.(place{3-cc}).(year_index)(exports.(place{cc}).(year_index)));
    firm_export_output.(place{cc}).(year_index) = [firm_export_output.(place{cc}).(year_index); ...
                            zeros(length(firm_output.(place{cc}).(year_index))-length(firm_export_output.(place{cc}).(year_index)),1)];    
    firm_export_revenue.(place{cc}).(year_index) = accumarray(varown_vec.(place{cc}).(year_index)(exports.(place{cc}).(year_index)),...
                            revenue.(place{3-cc}).(year_index)(exports.(place{cc}).(year_index))/import_tariff.(place{3-cc})); ...
    firm_export_revenue.(place{cc}).(year_index) = [firm_export_revenue.(place{cc}).(year_index); ...
                            zeros(length(firm_revenue.(place{cc}).(year_index))-length(firm_export_revenue.(place{cc}).(year_index)),1)];  
    firm_domestic_revenue.(place{cc}).(year_index) = firm_revenue.(place{cc}).(year_index)-firm_export_revenue.(place{cc}).(year_index);
    firm_domestic_output.(place{cc}).(year_index) = firm_output.(place{cc}).(year_index)-firm_export_output.(place{cc}).(year_index);
    firm_export_share.(place{cc}).(year_index) = firm_export_revenue.(place{cc}).(year_index)./firm_revenue.(place{cc}).(year_index);
    firm_TFPR.(place{cc}).(year_index) = firm_revenue.(place{cc}).(year_index)./firm_output.(place{cc}).(year_index);
    std_dev_firm_TFPR.(place{cc})(end+1) = nanstd(log(firm_TFPR.(place{cc}).(year_index)));
    var_firm_TFPR.(place{cc})(end+1) = nanvar(log(firm_TFPR.(place{cc}).(year_index)));
    
 
    % CATEGORIZE FIRMS INTO EXPORTERS AND NON-EXPORTERS
    exporters.(place{cc}).(year_index)=unique(varown_vec.(place{cc}).(year_index)(exports.(place{cc}).(year_index)));
    isnon_exporters.(place{cc}).(year_index)=~ismember(varown_vec.(place{cc}).(year_index),exporters.(place{cc}).(year_index))...
                                & ~isnan(varown_vec.(place{cc}).(year_index));
    non_exporters.(place{cc}).(year_index)=unique(varown_vec.(place{cc}).(year_index)(isnon_exporters.(place{cc}).(year_index)));
    num_non_exporters.(place{cc})=[num_non_exporters.(place{cc});...
                                   length(non_exporters.(place{cc}).(year_index))];    
    num_exporters.(place{cc})=[num_exporters.(place{cc});...
                                   length(exporters.(place{cc}).(year_index))];
    
    % NUMBER OF EXPORTED AND NONTRADED VARIETIES         
    num_exportED.(place{cc})=[num_exportED.(place{cc});...
                        sum(exports.(place{cc}).(year_index))];
    num_nontraded.(place{cc})=[num_nontraded.(place{cc});...
                    sum(nontraded.(place{cc}).(year_index))];
    
    % FIRM TFPR AND EMPLOYMENT MOMENTS          
    firm_TFPR_exporters.(place{cc}).(year_index) = firm_TFPR.(place{cc}).(year_index)(exporters.(place{cc}).(year_index));
    mean_firm_TFPR_exporters.(place{cc})(end+1) = nanmean(log(firm_TFPR_exporters.(place{cc}).(year_index)));
    std_dev_firm_TFPR_exporters.(place{cc})(end+1) = nanstd(log(firm_TFPR_exporters.(place{cc}).(year_index)));
    firm_TFPR_non_exporters.(place{cc}).(year_index) = firm_TFPR.(place{cc}).(year_index)(non_exporters.(place{cc}).(year_index));
    mean_firm_TFPR_non_exporters.(place{cc})(end+1) = nanmean(log(firm_TFPR_non_exporters.(place{cc}).(year_index)));
    std_dev_firm_TFPR_non_exporters.(place{cc})(end+1) = nanstd(log(firm_TFPR_non_exporters.(place{cc}).(year_index)));

    % MEAN EMPLOYMENT (NORMALIZED BY POPULATION)
    mean_employment_exporters.(place{cc})(end+1) = mean(firm_output.(place{cc}).(year_index)(exporters.(place{cc}).(year_index)) ...
        ./(sum(firm_output.(place{cc}).(year_index))));
    mean_employment_non_exporters.(place{cc})(end+1) = mean(firm_output.(place{cc}).(year_index)(non_exporters.(place{cc}).(year_index)) ...
        ./(sum(firm_output.(place{cc}).(year_index))));     
          
    % AVERAGE QUALITY STATS
    avg_qual_produced.(place{cc}).(year_index)=geomean(qualities_vec.(place{cc})(produced.(place{cc}).(year_index)));
    avg_qual_imported.(place{cc}).(year_index)=geomean(qualities_vec.(place{3-cc})(~produced.(place{cc}).(year_index)))*(omega^((-1)^(cc-1)));
	qual_consumed_vec.(place{cc}).(year_index) = [qualities_vec.(place{cc})(produced.(place{cc}).(year_index));...
                                                    qualities_vec.(place{3-cc})(~produced.(place{cc}).(year_index))];
    avg_qual_consumed.(place{cc}).(year_index)=geomean(qual_consumed_vec.(place{cc}).(year_index));                                                    
    if isnan(omega)
        avg_qual_consumed.(place{cc}).(year_index)=avg_qual_produced.(place{cc}).(year_index);
    end
    
     % Calculate group export statistics
%      group_export_share.(place{cc}).(year_index) = sum(reshape(exports.(place{cc}).(year_index), n_varieties/group_size{cc}, group_size{cc}))/sum(exports.(place{cc}).(year_index));
%      [leading_export_group_share.(place{cc}).(year_index), leading_export_group.(place{cc}).(year_index)] = max(group_export_share.(place{cc}).(year_index));
         
     
     
     
    
    % OTHER STATISTICS THAT CAN'T BE CALCULATED IN YEAR 1
     if year==1
        num_stopped_exporters.(place{cc})=NaN;
        num_new_exporters.(place{cc})=NaN;
        num_nontraded_lost_for.(place{cc})=NaN;
        num_nontraded_lost_dom.(place{cc})=NaN;
        num_exportED_lost_for.(place{cc})=NaN;
        num_exportED_to_nontraded.(place{cc}) = NaN;
        num_exportED_lost_dom.(place{cc})=NaN;
        num_produced_lost_for.(place{cc})=NaN;
        num_produced_lost_dom.(place{cc})=NaN;
        qual_consumed_growth_rate.(place{cc})=NaN;
        CES_growth_rate_vec.(place{cc}) = NaN;
        
        % VARIETY TRANSITIONS
         num_exported_to_exported.(place{cc})=NaN;
         num_exported_to_nontraded.(place{cc})=NaN;
         num_exported_to_imported.(place{cc})=NaN;
         num_nontraded_to_exported.(place{cc})=NaN;
         num_nontraded_to_nontraded.(place{cc})=NaN;
         num_nontraded_to_imported.(place{cc})=NaN;
         num_imported_to_exported.(place{cc})=NaN;
         num_imported_to_nontraded.(place{cc})=NaN;
         num_imported_to_imported.(place{cc})=NaN;
         
         prob_import_innovation_to_export.(place{cc})=NaN;
         prob_import_innovation_to_nontraded.(place{cc})=NaN;
         prob_nontraded_innovation_to_export.(place{cc})=NaN;
         prob_nontraded_innovation_to_nontraded.(place{cc})=NaN;
         
         CD_4_test.(place{cc})=NaN;
  else
         
         % GROWTH
         real_wage_growth_vec.(place{cc}) = [real_wage_growth_vec.(place{cc}); 
                        wage_level_vec.(place{cc})(end)/wage_level_vec.(place{cc})(end-1)];
         qual_consumed_growth_rate.(place{cc})=[qual_consumed_growth_rate.(place{cc});...
                                 avg_qual_consumed.(place{cc}).(year_index)/avg_qual_consumed.(place{cc}).(lag_year_index)];
    
         % MARKUPS STATS
         avg_markups_produced_vec.(place{cc}) = [avg_markups_produced_vec.(place{cc}); ...
             geomean([markups.(place{cc}).(year_index)(produced.(place{cc}).(year_index)); markups.(place{3-cc}).(year_index)(exports.(place{cc}).(year_index))])];
         avg_markups_consumed_vec.(place{cc}) = [avg_markups_consumed_vec.(place{cc}); ...
             geomean(markups.(place{cc}).(year_index))];
         dispersion_markups_consumed_vec.(place{cc}) = [dispersion_markups_consumed_vec.(place{cc}); ...
             std(log(markups.(place{cc}).(year_index)))];
         
         % EXPORTER TURNOVER - EXPORTERS THAT STOP EXPORTING
         num_stopped_exporters.(place{cc})=[num_stopped_exporters.(place{cc});...
                         sum(ismember(non_exporters.(place{cc}).(year_index),...
                                       exporters.(place{cc}).(lag_year_index)))];
         % EXPORTERS THAT STOP EXPORTING AND THAT EXIT ALTOGETHER
         num_all_stopped_exporters.(place{cc})=[num_all_stopped_exporters.(place{cc});...
                        sum(~ismember(exporters.(place{cc}).(lag_year_index),...
                                    exporters.(place{cc}).(year_index)))];
         all_stopped_exporters.(place{cc}).(year_index) = exporters.(place{cc}).(lag_year_index)(...
                        ~ismember(exporters.(place{cc}).(lag_year_index),...
                                    exporters.(place{cc}).(year_index)));
                                
         % INCUMBENTS THAT START EXPORTING
         num_new_exporters.(place{cc})=[num_new_exporters.(place{cc});...
                         sum(ismember(exporters.(place{cc}).(year_index),...
                                       non_exporters.(place{cc}).(lag_year_index)))];
         % ALL NEW EXPORTERS INCLUDING INCUMBENTS AND ENTRANTS
         all_new_exporters.(place{cc}).(year_index) = exporters.(place{cc}).(year_index)(~ismember(exporters.(place{cc}).(year_index),...
                                        exporters.(place{cc}).(lag_year_index)));
         num_all_new_exporters.(place{cc})=[num_all_new_exporters.(place{cc});...
                        sum(~ismember(exporters.(place{cc}).(year_index),...
                                        exporters.(place{cc}).(lag_year_index)))];
         
         
         
         % VARIETY TRANSITIONS - NONTRADED LOST TO THE FOREIGN COUNTRY                                           
         num_nontraded_lost_for.(place{cc})=[num_nontraded_lost_for.(place{cc});...
                                    sum(nontraded.(place{cc}).(lag_year_index) ...
                                                & ~produced.(place{cc}).(year_index))];
         % EXPORTS LOST TO THE FOREIGN COUNTRY
         num_exportED_lost_for.(place{cc})=[num_exportED_lost_for.(place{cc});...
                                    sum(exports.(place{cc}).(lag_year_index) ...
                                                & ~produced.(place{cc}).(year_index))];
         % EXPORTS THAT GO TO NONTRADED                                 
         num_exportED_to_nontraded.(place{cc}) = [num_exportED_to_nontraded.(place{cc});...
                                    sum(exports.(place{cc}).(lag_year_index) ...
                                                & nontraded.(place{cc}).(year_index))];
         % VARIETIES THAT SWITCH FROM NON-IMPORTS TO IMPORTS                                  
         num_produced_lost_for.(place{cc})=[num_produced_lost_for.(place{cc});...
                                    sum(produced.(place{cc}).(lag_year_index) ...
                                                & ~produced.(place{cc}).(year_index))];
         % NONTRADED LOST TO ANOTHER DOMESTIC FIRM
         num_nontraded_lost_dom.(place{cc})=[num_nontraded_lost_dom.(place{cc});...
             sum(produced.(place{cc}).(year_index) & nontraded.(place{cc}).(lag_year_index) & ... 
             was_stolen.(place{cc}))];
         % EXPORTS LOST TO ANOTHER DOMESTIC FIRM
         num_exportED_lost_dom.(place{cc})=[num_exportED_lost_dom.(place{cc}); ...
             sum(produced.(place{cc}).(year_index) & exports.(place{cc}).(lag_year_index) & ...
             was_stolen.(place{cc}))];
         % ALL VARIETIES STOLEN BY A DOMESTIC FIRM
         num_produced_lost_dom.(place{cc})=[num_produced_lost_dom.(place{cc});...
             sum(produced.(place{cc}).(year_index) & produced.(place{cc}).(lag_year_index) & ... 
             was_stolen.(place{cc}))];
         
         % WHAT HAPPENED TO EXPORTS?
         num_exported_to_exported.(place{cc}) = [num_exported_to_exported.(place{cc});...
                                sum(exports.(place{cc}).(year_index) & exports.(place{cc}).(lag_year_index))];
         num_exported_to_nontraded.(place{cc}) = [num_exported_to_nontraded.(place{cc});...
                                sum(nontraded.(place{cc}).(year_index) & exports.(place{cc}).(lag_year_index))];
         num_exported_to_imported.(place{cc}) = [num_exported_to_imported.(place{cc});...
                                sum(~produced.(place{cc}).(year_index) & exports.(place{cc}).(lag_year_index))];
         % WHAT HAPPENED TO NON-TRADED VARIETIES?
         num_nontraded_to_exported.(place{cc}) = [num_nontraded_to_exported.(place{cc});...
                                sum(exports.(place{cc}).(year_index) & nontraded.(place{cc}).(lag_year_index))];
         num_nontraded_to_nontraded.(place{cc}) = [num_nontraded_to_nontraded.(place{cc});...
                                sum(nontraded.(place{cc}).(year_index) & nontraded.(place{cc}).(lag_year_index))];
         num_nontraded_to_imported.(place{cc}) = [num_nontraded_to_imported.(place{cc});...
                                sum(~produced.(place{cc}).(year_index) & nontraded.(place{cc}).(lag_year_index))];
         % WHAT HAPPENED TO IMPORTED VARIETIES?
         num_imported_to_exported.(place{cc}) = [num_imported_to_exported.(place{cc});...
                                sum(exports.(place{cc}).(year_index) & ~produced.(place{cc}).(lag_year_index))];
         num_imported_to_nontraded.(place{cc}) = [num_imported_to_nontraded.(place{cc});...
                                sum(nontraded.(place{cc}).(year_index) & ~produced.(place{cc}).(lag_year_index))];
         num_imported_to_imported.(place{cc}) = [num_imported_to_imported.(place{cc});...
                                sum(~produced.(place{cc}).(year_index) & ~produced.(place{cc}).(lag_year_index))];
         
         % PROB OF EACH VARIETY CONDITIONAL ON INNOVATION
         prob_import_innovation_to_export.(place{cc}) = [prob_import_innovation_to_export.(place{cc}); ...
                                                                            sum(exports.(place{3-cc}).(lag_year_index) & was_stolen.(place{cc}) & exports.(place{cc}).(year_index)) ...
                                                                                                                ./ sum(exports.(place{3-cc}).(lag_year_index) & was_stolen.(place{cc}))];
         
         prob_import_innovation_to_nontraded.(place{cc}) = [prob_import_innovation_to_nontraded.(place{cc}); ...
                                                                            sum(exports.(place{3-cc}).(lag_year_index) & was_stolen.(place{cc}) & nontraded.(place{cc}).(year_index)) ...
                                                                                                                ./ sum(exports.(place{3-cc}).(lag_year_index) & was_stolen.(place{cc}))];
         
         prob_nontraded_innovation_to_export.(place{cc}) = [prob_nontraded_innovation_to_export.(place{cc}); ...
                                                                            sum(nontraded.(place{cc}).(lag_year_index) & was_stolen.(place{cc}) & exports.(place{cc}).(year_index)) ...
                                                                                                                ./ sum(nontraded.(place{cc}).(lag_year_index) & was_stolen.(place{cc}))];
         
         prob_nontraded_innovation_to_nontraded.(place{cc}) = [prob_nontraded_innovation_to_nontraded.(place{cc}); ...
                                                                            sum(nontraded.(place{cc}).(lag_year_index) & was_stolen.(place{cc}) & nontraded.(place{cc}).(year_index)) ...
                                                                                                                ./ sum(nontraded.(place{cc}).(lag_year_index) & was_stolen.(place{cc}))];
         
         CD_4_test.(place{cc}) = [CD_4_test.(place{cc}); ...
             sum(nontraded.(place{cc}).(year_index) & exports.(place{cc}).(lag_year_index) & was_stolen.(place{3-cc}))];
     end    
    
     
end